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Abstract 

We prove fixed points results for sandpiles starting with arbitrary ini- 
tial conditions. We give an effective algorithm for computing such fixed 
points, and we refine it in the particular case of SPM. 

1 Introduction 

Sandpiles are a simple but meaningful formal model for the simulation of systems 
governed by self-organized criticality (SOC). These systems, starting from any 
initial configuration, evolve to a "critical" state. Any perturbation of this critical 
state, no matter how small, originates a deep and uncontrollable reorganization 
of the whole. Then, they start evolving towards another critical state and so on. 
SOC systems are commonly used for simulating natural phenomena like snow 
avalanches, dune formations, but also woods fires and even, stock exchange 
crashes. 

The first discrete model for sandpiles, (later on) called SPM (Sand Pile 
Model), has been introduced in pQ. It is based on a simple local rule: a sand 
grain falls on its right if the difference between its sandpile and the one on its 
right is bigger than a certain amount of grains. In PJ|2|SJ|H], the model has 
been mathematically formalized and studied as a discrete dynamical system. In 
particular, they proved that SPM has fixed point dynamics and exhibited for- 
mulas for the precise expression of fixed points and for computing the transient 
length. 

Similar results exist for a more complete model, IPM(fc) (Ice Pile Model) 
introduced in The results found for these two models, synthesized in 0], 
are very interesting and complete but they concern only very special initial 
conditions in which all the sand grains are concentrated in a unique pile and 
there is no grain elsewhere. 

In this paper we generalize those results to arbitrary initial conditions. Of 
course, due to much greater combinatorial complexity of the problem, we do not 
have nice formulas but we give a fast algorithm for computing the fixed point. 

This paper is structured as follows. The next section recalls basic definitions 
about the main sandpile models. Section [21 resumes the main known results, 
which are generalized to arbitrary sandpiles in Section 0] Section [S] improves 
the results of the previous section for the special case of SPM. In the final section 
we draw our conclusion and give some perspectives. 
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2 Basic definitions 



A sandpile is a finite sequence of integers (ax, . . . , a;); £ S N is the length of 
the pile. Sometimes a sandpile is also called a configuration. Given a sandpile 
(oi, . . . , ai), the integer n = XL=i a» is the number of grains of the pile. Given a 
configuration (ai, . . . , a;), a subsequence Oj, . . . ,Oj (with 1 < z , j < Z and i < j) 
is a plateau if a*, = a^+i for i < k < j; & subsequence a^ai+i is a cZi/f if 
«i - flj+i > 2. 

In the sequel, each sandpile (a±, . . . , a;) will be conveniently represented on 
a two dimensional grid where Oj is the grain content of column i. 

A sandpile system is a finite set of rules that tell how the sandpile is updated. 
SPM is the most known and the most simple sandpile system. It consists in 
just one local rule. Moreover, all initial configurations contain n grains in the 
first column and nothing elsewhere i. e. they are of type (n) . The system rule 
can be defined in two equivalent ways. The former, introduced in yQ, considers 
the sequence Zi — ai — Oj+i of differences between two consecutive columns i 
and i + 1. If Zi > 2, then a sand grain falls from column i to i + 1 giving the 



following new sequence of differences (see Figure 1(a) I: 



z'i-i = Zi-i + 1 
z{ = z,-2 

4+i = z i+ i + i . 

The latter, introduced in 2^, deals with the real height of consecutive co- 
lumns. It has the advantage of having a simple and intuitive graphical repre- 



sentation. The updating rule is defined as follows (see Figure 1(b) I 

a ■ = a, - 1 
,/ _ „ ,i if Oi - a 4+ i > 2 . 

2 i+i — a i+i + 1 

Remark that there are also two ways of evolving the sandpile system: parallel 
and sequential execution mode. In parallel mode, all applicable rules are applied 
at once; only one rule at a time is applied in the sequential mode. For example, 
in Figure^] several grains might fall at the same time (columns 2 and 4 may lose 
a grain): the next configuration depends on the execution mode; if we choose 
the sequential mode, then the next configuration depends also on the column 
to which one applies the system rule. The execution mode is fixed from the 
beginning and does not change along the evolution of the system. In the sequel, 
we will be mainly interested in the sequential mode for the simplicity sake. 

IPM(fc) is another usual model that extends SPM. It contains the vertical 
rule of SPM defined above, plus a horizontal rule: the grains can slide on a 
horizontal plateau of length at most k as follows 

{...,p+l, p,p,...,p ,p-l,...) v (...,p,p,p,...,p,p,...) . 

V 

h' times (k' <k) 

This rule can be used to simulate slopes of less than 45° , while the vertical rule 
(when generalizing the definition of cliffs to any difference of height) simulates 
slopes bigger than 45°. 
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(a) Height differences. 




(b) Column grain content. 



Figure 1: updating in the SPM model. 



A sequence of configurations {ci} igN is called an orbit of initial condition ci 
if for all i G N, Cj+i is obtained from C{ via the application of a system rule. 
Remark that there might be more than one orbit for the same initial condition. 

The set of orbits O c with the same initial condition c is graphically repre- 
sented by the orbit graph Q c = (V,E), where V is the set of all configurations 
belonging to orbits in O c and (a, b) G E (E C V x V) if and only if b is obtained 
from a by an application of a system rule. Denote Q l c the orbit graph of the 
configuration c in which all configurations (including c) have length at most I. 
Given an orbit graph Q — (V,E) we say that a configuration c belongs to Q, 
denoted c G Q if c G V . A vertex v G V is a fixed point of a directed graph 
G = (V,E) if v has no outgoing edges. Given two graphs G\ = (Vi,Ei) and 
G2 = (V2, -B2)) Gi is a sub-graph of G 2 if Vi C V 2 and E 1 ! C E 2 - 

As an example, Figure El shows the orbit for the initial configuration (8) 
according to both sequential and parallel execution modes for SPM. One can 
remark that both sequential and parallel execution mode lead to the same fixed 
point (3,2,2,1). The difference between the two evolutions seems to consist 
only in the transient length. These remarks are true and will be justified in the 
following section. 
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(a) Sequential execution mode. 
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(b) Parallel execution mode. 



Figure 2: orbit of SPM for n = 8 with respect to execution mode. 

3 Known results 

The results of this section are essentially taken from |2j and (Hj. We stress that, 
for simplicity sake, only the sequential execution mode is considered. 



3.1 Fixed points 

First we recall the results for SPM, and for the slightly more complex model 
IPM(fc). 

The following theorem shows that, for any initial configuration (n), the orbit 
graph of SPM has a very special structure. 

Theorem 1 (j2]) For any integer n, the orbit graph Q{ n ) for SPM is a lattice 
and is finite. 

As a consequence of Theorem^we have that SPM has fixed point dynamics 
starting from any configuration (n). Moreover, this fixed point is unique and it 
is reached regardless of the order in which transitions are made and regardless of 
the execution mode (parallel or sequential) that has been chosen. The following 
lemma characterizes the elements of the lattice. 
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Lemma 2 (|_5J) Consider a configuration c and let n be its number of grains. 
Then, c € G(n) f or SPM if and only if it is decreasing and between any two 
plateaus there is at least a cliff. 

Remark 1 Consider a configuration c and let n be its number of grains. As- 
sume that c contains a plateau of length 3. Such a plateau can be seen as two 
consecutive plateaus of length 2. Thus, by Lemma\^ c does not belong to G( n )- 

The previous condition is not necessary for the characterization of the fixed 
point. Notwithstanding, it allows to obtain a much simpler proof of the following 
Theorem 

Notation. A couple of integers (p, k) is the decomposition of n <E N in its 
integer sum if n = k + i — k + Elfctli, 

Theorem 3 ( [2j ) There exists a unique decomposition of n 6 N in its integer 
sum. Then, the fixed point II obtained starting from the initial configuration (n) 
is the following 

n= f (p,p-l,...,l) ifk = 

\ (p,p — 1, . . . , k + 1, k, k, k — 1, . . . , 1) otherwise . 

Similar results hold for IPM(fc), we recall the main ones here. They can be 
found in their original form in [5]. 

Theorem 4 (|5j) For any integer n, the orbit graph Qua for IPM(k) is a lattice 
and is finite. 

Again, there exist a characterization of the configurations of the lattice |Hj • 
It is a generalization of Lemma|21which will not be used explicitly here. Remark 
that this allows to give the exact form of the fixed point of a configuration for 
IPM(fc). 

3.2 Transient length 

We have seen that according to the usual models, the sandpile (n) evolves to 
a fixed point. It can be interesting to know how much time {i.e. how many 
iterations of the system rule) it takes to reach such a fixed point - of course this 
time depends on the execution mode. 

For SPM, in the sequential execution mode, it is easy to compute the number 
of steps needed to reach the fixed point: it suffices to remark that the grains in 
the i-th column took i iterations to reach it; taking the sum all over the number 
of grains of the fixed point one finds the number of iterations. 

Theorem 5 ( j2] ) The length of the transient to reach the fixed point starting 
from the initial configuration (n) for SPM is given by the following formula 

t seq = \(P+ 1 )P(P ~ 1) + \ k ^V +l-k) = 0(n 3 / 2 ) , 
where (p, k) is the decomposition of n in its integer sum. 
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When parallel execution mode is used for SPM, things are a little bit more 
complex and one can only give an upper and a lower bound for the transient 
length. 

Theorem 6 ( [2J ) In the parallel execution mode, the length of the transient to 
reach the fixed point starting from the initial configuration (n) for SPM can be 
bounded as follows 

0(n) = < tpar < t seq = 0{n z ' 2 ) , 

where (p, k) is the decomposition of n in its integer sum. 

Finally for IPM(fc), all the paths from (n) to the fixed point do not neces- 
sarily have the same length, but the longest chain in sequential mode is known. 
Its exact expression can be found in 5 . 

4 Generalization to arbitrary sandpiles 

The results from the previous section are very complete, but they only apply to 
the particular case of initial configurations of length one. 

A little more general study was started in [2], where the authors remarked 
that when starting from decreasing configurations (for all < i < I, a t > a,; + i) 
the lattice structure is maintained with SPM. This result is extended in [H], 
where the author exhibits the set of fixed points but she does not associate to 
each initial configuration the corresponding fixed point. 

In this paper, we try to generalize these results to arbitrary initial configu- 
rations giving a fast algorithm for computing fixed points. 

Description of the algorithm. It consists in the iteration of two major steps: 

Cut: the configuration is divided into intervals so that the formulas of the next 
step can be applied; 

Compute: each interval is analyzed trying to figure out how it will be within a 
few iterations steps of the model. 

Remark that it is necessary to iterate since some grains can pass from 
one interval to another. Consider the case of two isolated columns of grains 
(m, 0, . . . , 0, n). The first column will collapse. Then, depending on the model 
chosen and on the value of m, n and of the gap between the two columns, the 
grains of the first column will be blocked by the second one or they will partially 
cover it. 

From now on, to simplify notations and proofs the results will be given only 
for SPM. Their extension to other models such as IPM(fc) is straightforward, 
using the formulas given in [oj. In fact our results extend to all models satisfying 
two conditions: lattice structure of the orbit graph, and reachability detection 
(one must be able to say whether a given configuration is reachable or not, in a 
way similar to Lemma|5J). Hence we will implicitly consider this class of sandpile 
models in the sequel of the paper when talking about "all" models. Next section 
justifies the restriction to these models. 
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4.1 Cut 



The way we are going to split the configuration into intervals is very simple: 
each interval has to contain a configuration which is reachable by the model, 
and it is the longest satisfying this property. 

At that point, the previous results need to be extended. Indeed, their authors 
supposed that the grains could move as far as possible. In the present case, the 
movement is limited to fixed intervals, grains are prevented from going too far 
on the right. All the results remain true, but have to be reformulated. 

Lemma 7 For all initial configurations c for SPM, Q c has no cycles. 

Proof. Given a configuration c = (oi, ...,a/), consider the quantity ip(c) = 

J2i=i Yl°jLi j- If c ' — ( a 'i; • ■ ■ 7 a [>) is obtained from c by applying the SPM rule 
at column i < I, then (p(c') = ip(c) — ai + a' i+1 = ip{c) — ai + a^+i + 1. Since the 
rule could be applied, ai > a; + i + 2, which implies ip(c') < y(c); in other words, 
the SPM system rule is irreversible. □ 

Our algorithm will only work with models having this irreversibility property. 
IPM(fc) has it, and any realistic model should have it. 

Lemma 8 For all initial configurations c of length at most I for SPM, Q l c has 
a unique fixed point II. Moreover, II 6 Q c . 

Proof. Let c be a configuration of length at most I, with n grains. First of all 
remark that Q l c is a sub-graph of Q c since for every configuration in Q l c , the path 
from the root (n) to this configuration is also in Q c (it just consists in applying 
system rules to columns of index less than I). 

Hence Q l c is finite, and by Lemma it has no cycles. Therefore every path 
starting from a configuration c reaches a fixed point. Assume that LIi and II2 
are two distinct fixed points of Q l c . They also belong to Q c , so they satisfy 
Lemma [21 i. e. they are decreasing and contain at most one plateau (as they 
are fixed points of SPM, there are no cliffs). Their structure is represented on 
Figure They are simply "staircases" with the addition of k grains: 

n f (p j ,p j -l,...,p j -l + l) if kj = 

3 I (P31P3 ,Pj ~ aj,P 3 -aj,...,pj-l + 2) otherwise, 

with Otj = I — kj ; — 1. 

Counting the number of grains in EL and II2 gives 

1 1 
n = ^2(pi + l-i)+h = ^2(P2 + 1 - i) + h (0 < h,k 2 < I) 

i=l i=l (1) 

= -l(2 Pl +l-t)+ki = ^l{2p 2 + l-l)+k 2 ■ 

This implies that |j>2 — Pi | = 7 • \k 2 — k\ | < 1, and as p\ and p 2 are integers we 
have pi = p 2 . It follows that fci = k 2 , and EL = U 2 . □ 

The irreversibility condition needed to apply the algorithm to a particular 
model can be relaxed. In fact, what is needed for the previous lemma to work 
is a lattice structure. This ensures irreversibility and unicity of the fixed point. 
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Figure 3: structure of a fixed point in Q c . The symbol q will be defined in 
Section Ol 

Proposition 9 For all initial configurations c of length at most I for SPM, Q l c 
is a lattice and Q l c C Q c . 

Proof. In the proof of Lemma |H1 we saw that Q l c is a sub-graph of Q c . By 
Lemma and |H1 we have the thesis. □ 

Lemma 10 Consider the SPM model, a configuration c of length at most I and 
let n be its number of grains. Then, c £ Gl n s if and only if it is decreasing and 
between any two plateaus there is at least a cliff. 

Proof. Consider a configuration c £ Gt n \ > where n is the number of grains of c. 
By Proposition |5J c belongs also to G( n ) an d hence it satisfies the hypothesis of 



For the opposite implication, let c — (ai, . . . , a^) be a decreasing configura- 
tion of length at most I in which any two plateaus have a cliff in between. Let n 
be the number of grain of c. By Lemma |2 c belongs to Gi n )- If we consider the 
path going up to the root, all the configurations encountered are of length less 
than I since the application of the system rule can only increase the length of a 
configuration. All the applications of the system rule take place in the interval 
[1, 1 — 1], hence each element of the path is also in Gr n y Q 

The previous lemma underlines the second condition that is needed for the 
sandpile model used in the algorithm: it should be possible to detect whether 
a configuration is reachable by the model. Moreover, for complexity issues (see 
Section POj) . it should be as fast as possible (for SPM and IPM(fc), it is linear, 
as a single scan suffices). 

Using the previous lemma one can complete the characterization of G l that 
begun with Proposition El 

Proposition 11 For any integer n it holds that 



where for any graph G = (V,E), G[V] is the sub-graph generated by the set of 
vertices V' C V. 



Lemma El 



G{ n) =G(n)[{ceN 1 }} 
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The previous proposition means that for SPM, . is exactly the sub-graph 
of G( n ) restricted to the configurations of length at most I, keeping all edges 
between these configurations. By Proposition this sub-graph is also a lat- 
tice. This holds for any model having the lattice structure and the reachability 
detection, this is the case for example for IPM(fc). 

The cut step is performed using Lemma Propositions 151 andll II A scan is 
performed in order to find the maximal intervals in which the corresponding sub- 
configurations are in Q c . For example, for SPM, a new interval starts whenever 
there are two consecutive columns i and i + 1 such that <Zj < Oj+i, or there is a 
second plateau not separated by a cliff (see Listing 



procedure cut (c[]) { // c is the initial configuration 
nbp =0; // number of plateaus 

1=0; // set of right extremities of the intervals 

for (i=l; i<l ; i++) { 

if (c[i+l] > c[i]) { // increase 

I = I U {i}; 
nbp = 0; 

} else if (c[i] - c[i + l] >= 2) { //cliff 
nbp = ; 

} else if (c[i + l] = c[i]) { // plateau 

nbp ++; 

if (nbp = 2) { 
I = I U {i}; 
nbp = ; 

} 

} 

} 

return I ; 

} 

Listing 1: procedure for cutting a configuration into intervals using SPM. 

By Lemma El the sub-configurations given by the intervals are in Q c and 
they are "maximal". Moreover, by Propositions 151 and 1111 we know that each 
sub-configuration is in G/ n .s and that it reaches the fixed point of G/ n .y where 

• Cj is the configuration corresponding to the i th interval, c; = (ak)kei ( > 

• li = | J, | is the length of c, ; 

• rii is the number of grains of Cj . 

Remark that the quantities U and rii can be computed inside the procedure cut. 

The last interval - reached when the scanning procedure arrives at a/ - is a 
"special case" : it is treated exactly like in the usual model, supposing there is 
as much space as necessary. An example of "cut" for the SPM model is shown 
in Figure 
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Figure 4: example of intervals after the cut step, for SPM. 



4.2 Compute 

In this section we show formulas for finding the fixed point locally to each 
interval. From now on and only for this section, we consider a generic interval 
and hence the index i will be omitted in the notations. Again, only SPM will be 
considered to simplify the formulas, but similar results hold for other models. 

The structure of the global configuration after the cut has to be precisely 
determined: we shall find the fixed point of every configuration in each interval. 

In the proof of Lemma |H1 we saw that the fixed point II of an interval with 
SPM is a kind of "staircase" as depicted in Figure El We need the expression of 
p (height of the leftmost column of the fixed point), q (height of the rightmost 
column) and k (number of exceeding grains) in terms of n and /. 

Equation implies that 



Then, k 



hl(2p + 1 — I) . For q, the same calculation is done, counting the 



number of grains in the fixed point 



^2(q + 1 - i) - k' (k' = l-k mod/, hence < k' < I) 



one finds q 



Remark 2 The quantity k is not necessarily a real number of grains. Indeed, 
when the fixed point finishes by ( . . . , 2 , 1 , , 0) , we have k = 1 . This does not cor- 
respond to a grain, but compensates a negative grain introduced when computing 
p (the last term of the sum equals —1). 

Moreover, the last interval has to be treated in another way, as it has infinite 
length. Choosing I — v' 8 "+ 1 ~ 1 , value obtained from the results in [2], would 
solve the problem. It corresponds to the length of the fixed point obtained from 
the configuration (n). With this new value, the calculation of p and q is correct. 
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4.3 Correctness 



In this section we will prove that the algorithm finishes and gives the correct 
fixed point. All results are stated for all suitable models. Proofs are given for 
the SPM case only, for simplicity sake. 

The superscript * for the quantities Ii,Ci,rii mean that we use their value 
they have at the end of the algorithm. 

Proposition 12 For any configuration, the algorithm finishes and returns a 
fixed point. 

Proof. Because of the irreversibility (Lemma EJ and of the finite number of 
grains, the algorithm finishes. 

Suppose that it does not return a fixed point, i.e. that it returns a config- 
uration c = (cij) where there is an index i such that a grain from et^ can move 
to a.i + \. Then ai > dj+i + 2, which implies that after the next cut step, Oj and 
Oi+i will belong to the same interval I (by definition of the cut). As there is 
a possible evolution in I, the algorithm has to compute the new configuration 
where at least this grain moved. □ 

We have to make sure that the fixed point characterized by the previous 
proposition is the "right" fixed point. 

The following lemma is quite straightforward for SPM. For other models, the 
proof depends on the model itself, although it comes from the irreversibility. For 
example for IPM(fc), there will be more cases to consider but the idea remains 
the same. 

Lemma 13 If a grain can move past a column i, then if any movement of 
grain passing this column is blocked, at least one of these movements will always 
remain possible. 

Proof. In the case of SPM the lemma can be reformulated as follows: a cliff at 
column k remains until the system rule is applied at column k. 

Suppose that the SPM system rule can be applied to a configuration c = 
(oi, . . . , a{) at column k. It means that a k > a k +\ + 2. If the rule applies at 
position j ^ k, and if we note d = (a'±, . . . , a' v ) the new configuration, there are 
four cases: 

• if j < k — 1, a' k = ak and a' k+l = dk+i; 

• i£j = k — l, a' k =a,k + l and a' k+l = ajt+i; 

• if j = k + 1, a' k = a k and a' k+1 = a k+ i - 1; 

• if j > k + 1, a' k = a k and a' k+l = a k+ i . 

In all cases, a' k > a k > a k +i + 2 > a' k+1 + 2. Hence the rule can still be applied 
at column k. □ 

Theorem 14 For any configuration c, Q c is a lattice and its fixed point coin- 
cides with the fixed point found by the algorithm. 
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Proof. Consider the intervals // obtained after the application of our algorithm 
to c . Split c according to these intervals (remark that we do not use the intervals 
given by the first cut step, we only use the final ones) into sub-configurations 
Ck = {a k ) keI f. 

No grain of any Cj ever moves to another Cj. Suppose it is not the case; we 
are going to simulate the behavior of the model, but inhibiting the applications 
of the system rule which move grains from an interval to another, raising a 
contradiction. Such a simulation ends when all the partial configurations Ck 
reach a fixed point. Let if be an interval in which Cj was once supposed to 
lose a grain. At this point of the simulation, the rule can still be applied in if 
(Lemma I13[) . Moreover, the fixed points of Cj and Cj+i correspond respectively 
to c{ and c{ +1 computed by the algorithm, because they belong to two different 
orbit graphs (lattices) and they do not receive nor lose any grain (we recall 
that the rules moving grains from one interval to another have been inhibited) . 
Therefore, the rule should also be applicable for our algorithm at the same 
column in //, which is not the case otherwise the algorithm would not have 
returned fLemma ll2fl . 

Therefore, the behavior of each Ck is not influenced by the others, and hence 
the set of reachable configurations can be obtained by joining all the possible 
behaviors for every jf . Thus Q c is a lattice since it is the product of all the 
i 1 

lattices Q* fs . 

(»i ) 

The fixed point found by the algorithm is necessarily a fixed point of Q c . 
Hence it coincides with the fixed point of Q c since Q c is a lattice. □ 

4.4 Complexity 

Our algorithm is a loop divided into two parts: the cut step and the compute 
step. In the cut step, the initial configuration is scanned and its complexity is 
Oil). Moreover, we can assume that the compute step is done in constant time 
for each interval, hence in 0(1) for the entire configuration (there are at most I 
intervals, for a strictly increasing configuration for example). 

The number of iterations of the loop is a little harder to evaluate. Intuitively, 
if there are many intervals, which means that the configuration is mostly non 
decreasing, the fixed point will be reached quite soon. If there are few intervals, 
then there will be few iterations because all the calculi will be done in the 
compute step. This is difficult to formalize in the general case, all we can do for 
now is give a quite large upper bound. 

Proposition 15 The algorithm performs at most \ -l- (l + 2f(n) — 1) iterations, 
where f(n) is the length of the fixed point reached by the configuration (n) (for 
example, f(ri) = \(y/8n + 1 - l)/2] = O(Vn) for SPM, see Remark^). 

Proof. First, note that there are at most I intervals. Consider the changes in 
the bounds of the intervals (and not in terms of grain content) between two 
iterations, i.e. between two cut steps. When there is no change, the algorithm 
returns. Hence at least one of the intervals "evolves" at each iteration, except 
the last one. 

Consider the upper bound (highest index) m of the interval Ii at the begin- 
ning of the program. Clearly, m > i, hence Ui can increase at most l + f(n) — i — l 
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times. It can increase one by one, until it reaches the end of the configuration 
or it disappears. I + f(n) — 1 is the maximum length of a configuration of size I 
with n grains, it is obtained when nearly all grains are in the last column. Since 
at least one Ui increases at each step (except the last step), there are at most 

l J+/(n)-2 

J2(l + f(n)-i-l)= i=2' l ' {l + 2f(n) ~ S) 

i=l i=/(n)-l 

changes of the intervals. 

Add I — 1 iterations for the loss of intervals, plus one final iteration to detect 
that nothing happened, and the result follows. □ 

Each iteration has a complexity of / (provided the cut step is linear and the 
compute step is constant, which is the case for SPM and IPM(fc)), hence the 
global complexity in the general case is in 0(l 2 ■ (I + 2f(n))). This is not so 
interesting, because it does not apply to any model in particular. Therefore it 
is not possible to give a sense to this result, comparing it to previous results. 

But it will be refined in the next section in the special case of SPM, and 
then it will be possible to compare it to the classical simulation. 

5 Improvement for SPM 

The algorithm described above works for the main existing sandpile models, 
provided the model has a lattice structure in which every element can be easily 
characterized. This ensures unicity of the fixed point, and the characterization 
enables to cut the configuration, to compute the fixed point faster. 
The simplicity of SPM allows a few optimizations. 

5.1 Merge 

The main optimization which can be achieved with SPM is the removal of the 
iteration over the cut step. The scan of the configuration can be replaced by a 
scan of the intervals in order to merge successive intervals when possible. The 
new structure of the algorithm is represented in Figure This operation does 
not affect the theoretical complexity, as the number of intervals can be as high 
as the number of columns, but in practice the gain is very important (very few 
intervals are actually of size 1, and they tend to disappear upon iterations). 

At each iteration of the merge step, the algorithm tries to combine intervals 
by considering the difference of height at their borders. If the merge succeeds, 
then new values are computed for the new interval. 

When looking at the border between two intervals Ii and h+i, there are two 
possibilities. 

• Either qi < Pi+\ + 1, in which case nothing can happen at the border. 
There are no cliffs inside the intervals by construction of the fixed point, 
hence nothing has to be done. 

• Or qi > Pi+\ + 2, which means that the system rule can be applied at the 
border. To obtain a new fixed point, the two intervals will be merged into 
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Figure 5: structure of the algorithm refined for SPM 

I' a (j\ = I-i U Ii+i, where a is some suitable renumbering function. Remark 
that from one iteration to another, some intervals disappear (merge) and 
some remain, introducing lag in the indices. 

By Lemma lTUl in this new interval the configuration belongs to @u^+n* t)' 
Indeed, there is at most one plateau in Jj, at most another one in 
and they are separated by the cliff at the border. Therefore, the previous 
formulas still hold in the new interval I a U) with 

rii + n l+ i 
h + U+i 

i> 

Due to the fact that we are in the lattice GiZi \, we know that the com- 

Q(i .' 

puted fixed point will be correct. Therefore this operation can efficiently 
merge two intervals. 

Once two intervals have merged, the interval list has to be updated with the 
new interval which replaces the two old ones. That way, the next iteration will 
act on the new list. This allows a newly created interval to merge at the next 
step, letting grains move as far as possible. 

5.2 Transient length 

For SPM, the algorithm can also compute the transient length. It is not possible 
in general because all paths from the root to the fixed point may not have the 
same length (IPM(fc) for example), but it is the case with SPM. So the number 
of steps simulated by the algorithm will correspond to the transient length, 
independently from the choice of the intervals. 

To compute its value inside an interval, we associate to every sand grain a 
number corresponding to its movement, as in [2]. All these elementary values 



'a(i) 

l a(i) 

C 'aU) 
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are added, supposing the configuration in Figure [21 has been reached. Then we 
have to subtract the "movement" weight t of the initial configuration. One 
finds 

t=5>-o+ E i- f o = ^- 1 )^- 2f + 1 ) + ^- fc - 1 )- f o , ( 2 ) 

«= i=l— k 

where t° = X)'=o * ' a *+i can computed during the cut step; j is the index of 
the leftmost column of the interval under consideration. 

Remark 3 Again, the quantity k is not really a number of grains. The right- 
most sum of Equation can be too big, but this will be compensated by a 
smaller value on the left. The extra grain added on the right is counted as a 
negative grain on the left. 

One has to add all these values, over each interval at each iteration, to 
obtain the total transient length. When intervals do not merge, t does not 
change. Otherwise, we are going to count the number of steps simulated, and 
subtract everything that had already been done during last iteration. 

t' r , = -V ... (l' (V ... - 21' + l) +-k' (21' -k' -l) 

— U — — ■ h 

The last term is the number of steps spent to move Uj+x grains from interval Ii 
to 

5.3 Summary 

Table^sums up all the calculations that have to be computed at each iteration, 
in the SPM example. The notation f3j represents the value of the variable /3 in 
the i th interval at the j th iteration. 

At the end of the execution, these values allow to determine exactly the 
shape of the fixed point (p, q and k in each interval), and the total transient 
length. 

5.4 Complexity (for SPM) 

As seen above, the theoretical complexity is not improved with the modified 
algorithm, but we can refine the computation made in Section 14.41 

The compute and merge steps both consist in a scan of the intervals, hence 
they are in 0(1). The number of iterations is clearly bounded by 0(1), because 
there can not be more merging than the number of intervals. Moreover, the 
following proposition shows that it is also bounded by the number of grains. 
This last constraint becomes interesting when grains are scattered along the 
configuration, i.e. when / ^> n. 

Proposition 16 Consider a configuration with n grains. The algorithm finds 
the fixed point performing at most 0{n) merge steps. 
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Table 1: summary of the operations performed by the algorithm modified for 
SPM, during one iteration. Column A contains the values obtained when q\ < 
pl +1 + 1, and B the values for q\ > p 3 i+1 + 2. 

Proof. It suffices to prove that if there are m mergings, then there are at least 
m grains in the configuration. Indeed, in order for /; and Ji+i to merge we 
shall have qi > Pi+\ + 2 > 2. This means that It has at least two grains in its 
rightmost column, mark the lowest one. After the merging, the marked grain 
does not move because it is at height 1; moreover, it is not anymore at the border 
of two intervals. Therefore, any successive merging will mark a different grain, 
which has not been marked yet. Thus, there are no more than n mergings. □ 

The following example shows that the bound given in Proposition^] can be 
reached. 

Example 1 Consider a configuration c of length I defined as follows 

( 7 ifi = 4j , 0< j<A[n/7\ 

Vi € {1, ...,/}, Ci = < n mod 7 ifi = 4[n/7\ 

[ otherwise. 

In other words, c — (7, 0, 0, 0, 7, 0, 0, 0, . . . , 7, 0, 0, 0, n mod 7) (figure^). The 
cut step produces 2\n/7~\ — 1 intervals of type (7, 0, 0) and (0). After the compute 
step, the configuration is c' = (3, 2, 2, 0, 3, 2, 2, 0, . . . , 3, 2, 2, 0, X\, X2, Xa). One 
can easily see that there will be \n/7~\ — 1 mergings of intervals of type (3,2,2) 
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and (0) (in red, dashed dotted lines), in order to obtain the fixed point (3, 2, 1, 1, 
3, 2, 1,1,..., 3, 2,1,1, xi, x 2 ,x 3 ). 



Figure 6: example of configuration reaching the 0(n) complexity 

Our algorithm computes the fixed point of any initial configuration in time 
0(1 ■ min(n, I)). There is a gain of at least y/n compared to the naive simula- 
tion, whose complexity is 0(1 ■ n 3 / 2 ) (n 3 / 2 is the number of avalanches, while I 
corresponds to the scan of the configuration in order to find a cliff). 



6 Conclusions 

In this paper we proposed an algorithm for finding the fixed points of a large 
class of sandpile models starting from arbitrary finite configurations. 

Time complexity of the algorithm depends on the structure of the initial 
configuration and, of course, on the model chosen. It might be interesting 
to build hierarchies of models classified according to time complexity of the 
(specialized version of the) algorithm, and study the algebraic property of the 
induced orbit graphs. 

Another research direction consists in the study of models which allow grains 
to move in more than one direction and in parallel execution mode. Remark 
in fact that grains in a sandpile are essentially subject to two different type 
of forces: a vertical one due to gravity, and a horizontal one due to wind. 
Realistic models should consider both forces, and apply them in parallel to 
every grain. Our algorithm would be adapted to such models, and would allow 
much faster computation of this kind of natural phenomenon. Our research 
currently consists in finding such a model, with the needed lattice structure. 
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